%To generate CM4 predicted geomagnetic fields at observatories
% from 1995 to 2006. As part of the NASA project (See Gary Egbert's email
% on June 2, 2009.


%first load the adjusted solar flux data (spidr)
load c:\manoj\projects\gary\f107\Adjusted_Solar_Flux_f10.7_spidr_1995-2006.mat;
fday_solar_flux = datenum(Time_Date_stamp(:,1)) - datenum(2000,1,1);

%load the observatory data
load C:\Manoj\projects\gary\observatory_data_secular_removed.mat phi_gg station_name fday;;
load c:\manoj\projects\gary\corrected_theta_gg  theta_gg;
dummy =  repmat(f107_1995_2006, [1,24]);
solar_flux_hourly_new =  reshape(dummy', [length(f107_1995_2006)*24,1]);
%dst
load c:\manoj\projects\gary\dst\Dst_XXXX.txt;

fday_dst_all = Dst_XXXX(:,1);

start_fday = datenum(1995,1,1) - datenum(2000,1,1);
end_fday = datenum(2006,12,31,23,59,0) - datenum(2000,1,1);

L = fday_dst_all >= start_fday & fday_dst_all <= end_fday;

Dst_selected = Dst_XXXX(L,2:25);
fday_dst_selected = fday_dst_all(L);

dst_hourly = reshape(Dst_selected',[length(Dst_selected)*24,1]);
[a,b] = xlsread('c:\manoj\projects\gary\Gmagobs_heights.xls');
[Y,IA,IB] =  intersect(station_name,b(:,2));
height = a(IB,end);
L = isnan(height);
height(L) = 0; %heights verified against that at BGS geomag website
height_in_km = height./1000;%km



%run Cm4

%       path(1)='../umdl.CM4'       ! path for model coefficients
%       unit(1)=11                  ! logical unit for model coefficients
%       load(1)=.true.              ! set to initially read model coefficients
%       indx(1)=.false.             ! get Dst from argument
%       indx(2)=.false.             ! get solar flux from argument
%       gmut   =.true.              ! compute MUT from UT
%       cord   =.false.             ! use geocentric coordinates
%       pred(1)=.true.              ! compute internal B fields
%       pred(2)=.true.              ! compute magnetospheric B field
%       pred(3)=.true.              ! compute ionospheric B field
%       pred(4)=.true.              ! compute toroidal B field for...
%       pred(5)=.true.              ! ...Magsat
%       pred(6)=.true.              ! ...dawn LT
%       curr   =.false.             ! no current density J fields computed
%       coef   =.false.             ! no coefficients generated
%       nhmf(1)=13                  ! n_max for "core field"
%       nhmf(2)=40                  ! n_max for "crustal field"
%       nlmf(1)=1                   ! n_min for "core_field"
%       nlmf(2)=14                  ! n_min for "crustal field"
%       ut     =1980.d0             ! constant UT (years)
%       phi    =0.d0                ! constant longitude
%       alt    =430.d0              ! constant altitude
%       dst    =-10.d0              ! constant Dst
%       f107   =200.d0              ! constant solar flux
%       perr   =1                   ! error print flag
%       oerr   =6                   ! error logical unit
%       cerr   =0                   ! error flag




a = 6371.2;
t = fday + (30/(60*24)) ; % magnetic data and dst values are value over (i-1)th hour, thus values are
%centered on the half-hour. Currently fday(1) is 1995 1 1 0 0 0, where are it
%should be 1995 1 1 0 30 0 . +(30/(60*24))  will make this difference 
% note the results made on June 3-5, 2009 were not corrected for this

Dst = dst_hourly;
F107 = solar_flux_hourly_new;

options.n_min = [ 1 14]; % n_min
options.n_max = [13 45]; % n_max

% options.n_min = [1 1]; % n_min
% options.n_max = [0 0]; % n_max

options.cord = 1;     % geocentric coordinates
options.pred = [1 1 1 1 1 1];
%options.pred = [0 1 1 0 0 0]; %compute only magntospheric and ionospheric fields
options.filename =  'umdl.CM4';

%Compute for station NGK

for i = 29:160,

theta = theta_gg(i);
phi = phi_gg(i);
r = 6371.2 + height_in_km(i);

[B_x, B_y, B_z] = cm4(t, r, theta, phi, Dst, F107, options);

eval(['save ' 'c:\manoj\projects\gary\results\' sprintf('%d_%s',i,station_name(i,:))]);
display(['save ' 'c:\manoj\projects\gary\results\' sprintf('%d_%s',i,station_name(i,:))]);
end;



%display the results
load C:\Manoj\projects\gary\observatory_data_secular_removed.mat X_secular_removed Y_secular_removed Z_secular_removed station_name fday;
%%

%L = fday+datenum(2000,1,1) > datenum(2000,8,15) & fday+datenum(2000,1,1) < datenum(2000,9,15);
L = boolean(ones(size(L)));
for i = 1:160,
eval(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:)) ' B_x B_y B_z']);
display(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_new_%sB_x',i,station_name(i,:))]);
subplot(311);
plot(fday(L)+datenum(2000,1,1),X_secular_removed(L,i));
hold on;
%plot(fday(L)+datenum(2000,1,1),sum(B_x(L,3:7),2)-nanmean(sum(B_x(L,3:7),2)),'r');
plot(fday(L)+datenum(2000,1,1),sum(B_x(L,3:7),2),'r');
hold off;
title(['1995 - 2006' station_name(i,:)]);
axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);
datetick('x','keeplimits');
subplot(312);
plot(fday(L)+datenum(2000,1,1),Y_secular_removed(L,i));
hold on;
%plot(fday(L)+datenum(2000,1,1),sum(B_y(L,3:7),2)-nanmean(sum(B_y(L,3:7),2)),'r');
plot(fday(L)+datenum(2000,1,1),sum(B_y(L,3:7),2),'r');
axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);
hold off;
datetick('x','keeplimits');
subplot(313);
plot(fday(L)+datenum(2000,1,1),Z_secular_removed(L,i));
hold on;
%plot(fday(L)+datenum(2000,1,1),sum(B_z(L,3:7),2)-nanmean(sum(B_z(L,3:7),2)),'r');
plot(fday(L)+datenum(2000,1,1),sum(B_z(L,3:7),2),'r');
axis([datenum(1995,1,1) datenum(2006,12,31) -inf inf]);
hold off;

datetick('x','keeplimits');
%pause;
saveas(gcf,['C:\Manoj\projects\gary\cm4figures\' station_name(i,:)],'pdf');
end;


%%
load C:\Manoj\projects\gary\observatory_data_secular_removed.mat  station_name ;

B_x_Core = zeros([105192,160]);
B_x_Crustal = zeros([1,160]);
for i = 1:160,
eval(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:)) ' B_x']);
display(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:))]);
B_x_Core(:,i) = B_x(:,1);
B_x_Crustal(i) = B_x(1,2);
B_x_Dst_Primary(:,i) = B_x(:,3);
B_x_Dst_Secondary(:,i) = B_x(:,4);
B_x_Iono_Primary(:,i) = B_x(:,5);
B_x_Iono_Secondary(:,i) = B_x(:,6);
B_x_Toroidal(:,i)  = B_x(:,7);
clear B_x;
end;
save c:\manoj\projects\gary\B_x_Dst_Primary B_x_Dst_Primary;
save c:\manoj\projects\gary\B_x_Dst_Secondary B_x_Dst_Secondary;
save c:\manoj\projects\gary\B_x_Iono_Primary B_x_Iono_Primary;
save c:\manoj\projects\gary\B_x_Iono_Secondary B_x_Iono_Secondary;
save c:\manoj\projects\gary\B_x_Iono_Toroidal B_x_Toroidal;
save c:\manoj\projects\gary\B_x_Core B_x_Core
save c:\manoj\projects\gary\B_x_Crustal B_x_Crustal
%%
load C:\Manoj\projects\gary\observatory_data_secular_removed.mat  station_name ;
B_y_Core = zeros([105192,160]);
B_y_Crustal = zeros([1,160]);
B_y_Dst_Primary = B_y_Core;
B_y_Dst_Secondary = B_y_Core;
B_y_Iono_Primary = B_y_Core;
B_y_Iono_Secondary = B_y_Core;
B_y_Toroidal = B_y_Core;


for i = 1:160,
eval(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:)) ' B_y']);
display(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:))]);
B_y_Core(:,i) = B_y(:,1);
B_y_Crustal(i) = B_y(1,2);
B_y_Dst_Primary(:,i) = B_y(:,3);
B_y_Dst_Secondary(:,i) = B_y(:,4);
B_y_Iono_Primary(:,i) = B_y(:,5);
B_y_Iono_Secondary(:,i) = B_y(:,6);
B_y_Toroidal(:,i)  = B_y(:,7);
clear B_y;
end;

save c:\manoj\projects\gary\B_y_Dst_Primary B_y_Dst_Primary;
save c:\manoj\projects\gary\B_y_Dst_Secondary B_y_Dst_Secondary;
save c:\manoj\projects\gary\B_y_Iono_Primary B_y_Iono_Primary;
save c:\manoj\projects\gary\B_y_Iono_Secondary B_y_Iono_Secondary;
save c:\manoj\projects\gary\B_y_Toroidal B_y_Toroidal;
save c:\manoj\projects\gary\B_y_Core B_y_Core;
save c:\manoj\projects\gary\B_y_Crustal B_y_Crustal;

%%
load C:\Manoj\projects\gary\observatory_data_secular_removed.mat  station_name ;
B_z_Core = zeros([105192,160]);
B_z_Crustal = zeros([1,160]);
B_z_Dst_Primary = B_z_Core;
B_z_Dst_Secondary = B_z_Core;
B_z_Iono_Primary = B_z_Core;
B_z_Iono_Secondary = B_z_Core;
B_z_Toroidal = B_z_Core;


for i = 1:160,
eval(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:)) ' B_z']);
display(['load ' sprintf('C:\\Manoj\\projects\\gary\\results\\%d_%s',i,station_name(i,:))]);
B_z_Core(:,i) = B_z(:,1);
B_z_Crustal(i) = B_z(1,2);
B_z_Dst_Primary(:,i) = B_z(:,3);
B_z_Dst_Secondary(:,i) = B_z(:,4);
B_z_Iono_Primary(:,i) = B_z(:,5);
B_z_Iono_Secondary(:,i) = B_z(:,6);
B_z_Toroidal(:,i)  = B_z(:,7);
clear B_z;
end;

save c:\manoj\projects\gary\B_z_Dst_Primary B_z_Dst_Primary;
save c:\manoj\projects\gary\B_z_Dst_Secondary B_z_Dst_Secondary;
save c:\manoj\projects\gary\B_z_Iono_Primary B_z_Iono_Primary;
save c:\manoj\projects\gary\B_z_Iono_Secondary B_z_Iono_Secondary;
save c:\manoj\projects\gary\B_z_Iono_Toroidal B_z_Toroidal;
save c:\manoj\projects\gary\B_z_Core B_z_Core;
save c:\manoj\projects\gary\B_z_Crustal B_z_Crustal;
